home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / map_image.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  12.0 KB  |  309 lines
  1. ; $Id: map_image.pro,v 1.16 1997/03/17 23:40:06 dave Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. function map_image_missing, image_orig, max_value, min_value
  7. ;
  8. ; Return an array of 1's where the data are outside the range of min_value
  9. ; to max_value.  Max_value and/or min_value may be undefined.  If both are
  10. ; undefined, return a -1.
  11. ;
  12. k = n_elements(max_value) * 2 + n_elements(min_value)
  13. case k of
  14. 0: return, -1
  15. 1: return, Image_orig le min_value
  16. 2: return, Image_orig ge max_value
  17. 3: return, (Image_orig ge max_value) or (Image_orig le min_value)
  18. endcase
  19. end
  20.  
  21.  
  22.  
  23. FUNCTION  map_image, Image_Orig, Startx, Starty, xsize, ysize, $
  24.                 LATMIN = latmin,   LATMAX = latmax,  $
  25.                 LONMIN = lonmin,   LONMAX = lonmax,  $
  26.                 COMPRESS = compress, BILINEAR = bilin, $
  27.         WHOLE_MAP = whole_map, SCALE = scalef, $
  28.         MISSING = missing, MAX_VALUE = max_value, MIN_VALUE=min_value
  29. ;+NODOCUMENT
  30. ;NAME:
  31. ;     map_image
  32. ;PURPOSE:
  33. ;       This function returns the Image_Orig image warped to fit
  34. ;       the current map. Image_Orig must be centered at 0.  This
  35. ;     routine works in image space.
  36. ;Category:
  37. ;        Mapping
  38. ;Calling Sequence:
  39. ;        result = map_image(Image_Orig [, Startx, Starty [, xsize, ysize]])
  40. ;INPUT:
  41. ;      Image_Orig--- A two-dimensional array representing geographical
  42. ;               image to be overlayed on map.  It has Nx columns,
  43. ;        and Ny rows.
  44. ;KEYWORDS:
  45. ;    LATMIN --- the latitude corresponding to the first row of Image_Orig.
  46. ;        The default value is -90.  Latitude and Longitude values
  47. ;        refer to the CENTER of each cell.
  48. ;    LATMAX --- the latitude corresponding to last row of Image_Orig. The
  49. ;        default is  90 - (180. / Ny).
  50. ;    LONMIN --- the longitude corresponding to the first column of
  51. ;        Image_Orig. The default value is -180.
  52. ;    LONMAX --- the longitude corresponding to the last column
  53. ;        of Image_Orig. The default is  180. - (360./Nx).
  54. ;        For images crossing the international dateline, lonmax
  55. ;        will be less than lonmin.
  56. ;        If the longitude of the last column is equal to 
  57. ;        (lonmin - (360. /Nx)) MODULO 360.
  58. ;        it is assumed that the image covers all longitudes.
  59. ;    BILINEAR --- A flag, if set, to request bilinear interpolation. The
  60. ;        default is nearest neighbor.  Bilinear appears much better.
  61. ;    COMPRESS --- Interpolation compression flag.  Setting this to
  62. ;        a higher number saves time --- lower numbers produce
  63. ;        more accurate results.  Setting this to 1
  64. ;        solves the inverse map transformation for every
  65. ;        pixel of the output image.  Default = 4 for output devices
  66. ;        with fixed pixel sizes. Fix is used to make this an int.
  67. ;    SCALE = pixel / graphics scale factor for devices with scalable
  68. ;        pixels (e.g. PostScript).  Default = 0.02 pixels/graphic_coord.
  69. ;        This yields an approximate output image size of 350 x 250.
  70. ;        Make this number larger for more resolution (and larger
  71. ;        PostScript files, and images), or smaller for faster
  72. ;        and smaller, less accurate images.
  73. ;    MISSING = value to set areas outside the valid map coordinates.
  74. ;        If omitted, areas outside the map are set to 255 (white) if
  75. ;        the current graphics device is PostScript, or 0 otherwise.
  76. ;    MAX_VALUE = values in Image_Orig greater than or equal to MAX_VALUE
  77. ;        are considered missing.  Pixels in the output image
  78. ;        that depend upon missing pixels will be set to MISSING.
  79. ;    MIN_VALUE = values in Image_Orig less than or equal to MIN_VALUE
  80. ;        are considered missing.
  81. ; Optional Output Parameters:
  82. ;    Startx --- the  x coordinate where the left edge of the image
  83. ;        should be placed on the screen.
  84. ;    Starty --- the y coordinate where th bottom edge of the image
  85. ;        should be placed on the screen.
  86. ;    xsize ---  returns the width of the resulting image expressed in
  87. ;        graphic coordinate units.  If current graphics device has
  88. ;        scalable pixels,  the value of XSIZE and YSIZE should
  89. ;        be passed to the TV procedure.
  90. ;    ysize ---  returns the pixel height of the resulting image, if the
  91. ;        current graphics device has scalable pixels. 
  92. ;
  93. ;Output:
  94. ;      The warped image is returned.
  95. ;
  96. ; Procedure:  An image space algorithm is used, so the time required
  97. ;    is roughly proportional to the size of the final image.
  98. ;    For each pixel in the box enclosed by the current window,
  99. ;    and enclosed by the Image, the inverse coordinate transform is
  100. ;    applied to obtain lat/lon.  The lat/lon coordinates are then scaled
  101. ;    into image pixel coordinates, and these coordinates are then
  102. ;    interpolated from Image values.
  103. ;
  104. ; Restrictions: Probably won't work on VMS cause it doesn't support
  105. ;     NaN.
  106. ;
  107.  
  108. ;MODIFICATION HISTORY:
  109. ;       CAB, Feb, 1992. Map_image has been changed to handle images
  110. ;           crossing the international dateline in a more convenient way.
  111. ;           Specifically, it no longer requires that the keyword LONMIN be
  112. ;           greater than or equal to -180 or the keyword LONMAX be
  113. ;        less than or equal to 180.
  114. ;    DMS, Aug, 1992.  Completly rewritten.  Uses different algorithms.
  115. ;    DMS, Dec, 1992.  Coordinates were off by part of a pixel bin.
  116. ;        Also, round when not doing bi-linear interpolation.
  117. ;    DMS, Sep, 1993.  Added MAX_VALUE keyword.
  118. ;    DMS, Nov, 1994.  Added MIN_VALUE keyword.
  119. ;       SVP, Mar, 1995.  Compress is now fix()'d. y is now scaled correctly.
  120. ;       SVP, May, 1996.  Removed Whole_Map keyword. Changes in the noborder
  121. ;                        behavior of MAP_SET make this keyword obselete.
  122. ;    DMS, Nov, 1996.    Adapted for new maps, rev 2.
  123. ;-
  124.  
  125. ON_ERROR,2
  126.  
  127. ;t0 = systime(1)
  128. if (!x.type NE 2) and (!x.type ne 3) THEN  $        ;Need Mapping Coordinates
  129.    message, "Current window must have map coordinates"
  130.  
  131. s = size(Image_Orig)
  132. if s[0] ne 2 THEN message, " Image must be a two- dimensional array."
  133. s1 = s[1]           ; # of columns
  134. s2 = s[2]           ; # of rows
  135. if s[1] le 1 or s[2] le 1 THEN $
  136.     message, 'Each dimension must be greater than 1."
  137.  
  138. ; If both latmin & latmax are missing, assume image covers entire globe.
  139. if N_ELEMENTS(latmin) eq 0 then latmin = -90.
  140. if N_ELEMENTS(latmax) EQ 0 THEN latmax = 90. ;This is as documented
  141. ; if N_ELEMENTS(latmax) EQ 0 THEN latmax = 90. - 180./s2  ;This isn't
  142.  
  143. ; If both lonmin & lonmax are missing, assume image covers all longitudes
  144. ;       with duplication.
  145. if N_ELEMENTS(lonmin) EQ 0 THEN lonmin = -180.
  146. if N_ELEMENTS(lonmax) EQ 0 THEN lonmax = 180. ;This is as documented
  147. ; if N_ELEMENTS(lonmax) EQ 0 THEN lonmax = 180.- (360./s1)   ;This isn't
  148.  
  149. latmin = float(latmin)     &    lonmin = float(lonmin)
  150. latmax = float(latmax)     &    lonmax = float(lonmax)
  151.  
  152. ;    Does image wrap around globe?
  153. wrap = ((lonmin - 360./s1 + 720.) mod 360.) - ((lonmax + 720.) mod 360.)
  154. wrap = abs(wrap) lt 1e-4    ;Allow for roundoff
  155. ltlim = [latmin, latmax]
  156.  
  157. if wrap eq 0 then lnlim = [lonmin, lonmax] else lnlim = [-180., 180.]
  158.  
  159. ; Find the extent of the our limits in the map on the screen by
  160. ;       making a n x n array of lon/lats spaced over the extent of
  161. ;       the image and saving the extrema.
  162.  
  163.  
  164. scale = !d.flags and 1        ;TRUE if device has scalable pixels
  165. if n_elements(scalef) le 0 then scalef = 0.02   ;PostScript scale factor
  166. IF scale THEN BEGIN        ;Fudge for postscript?
  167.     !x.s = !x.s * scalef
  168.     !y.s = !y.s * scalef
  169. ENDIF ELSE scalef = 1
  170.  
  171. if n_elements(compress) le 0 then compress = 4   $ ;Default value
  172. else compress = fix(compress)
  173.  
  174. ; Missing data value should equal the background or user-supplied value.
  175. if n_elements(missing) le 0 then begin
  176.   if (!d.flags and 512 ne 0) then missing = !d.n_colors-1 else missing = 0
  177.   endif
  178.  
  179. screen_x = long(scalef * !x.window * !d.x_size) ;Map extent on screen
  180. screen_y = long(scalef * !y.window * !d.y_size)
  181.  
  182. ; See if we can use a smaller area than the plot window
  183. if wrap eq 0 and abs(latmax-latmin) gt 90 then begin
  184.     n = 31                      ;Subdivisions across lat/lon range.
  185.     x = (findgen(n) * ((lnlim[1]-lnlim[0])/float(n-1)) + lnlim[0]) # $
  186.       replicate(1.,n)
  187.     y = replicate(1.,n) # $
  188.       (findgen(n) * ((ltlim[1]-ltlim[0])/float(n-1)) + ltlim[0])
  189.     x = convert_coord(x,y, /DATA, /TO_DEVICE) ;Latlon to device
  190.     
  191.     y = reform(x[1,*], n^2, /OVER) ;Get device coords separately
  192.     x = reform(x[0,*], n^2, /OVER)
  193.     good = where(finite(x))
  194.     x = x[good] & y = y[good]
  195.     screen_x = long([screen_x[0] > min(x, MAX=maxx), screen_x[1] < maxx])
  196.     screen_y = long([screen_y[0] > min(y, MAX=maxy), screen_y[1] < maxy])
  197. endif
  198.  
  199. ;       Get next larger multiple of compress for resulting image size.
  200. nx = (screen_x[1] - screen_x[0]+compress)/compress
  201. ny = (screen_y[1] - screen_y[0]+compress)/compress
  202.  
  203.                         ;X and Y screen coordinates
  204. x = (findgen(nx) * compress) # replicate(1.0,ny) + screen_x[0]
  205. y = replicate(1.0, nx) # (findgen(ny) * compress) + screen_y[0]
  206.  
  207. x = convert_coord(x, y, /DEVICE, /TO_DATA)      ;Screen to lat/lon
  208. y = reform(x[1,*], nx, ny, /OVER)       ;Separate lat/longit
  209. x = reform(x[0,*], nx, ny, /OVER)
  210.  
  211. bad = where(finite(x) eq 0b, count) ;Not all machines handle NaN properly
  212. if count ne 0 then begin        ;So fake it for points off the map.
  213.     x[bad] = 1.0e10
  214.     y[bad] = 1.0e10
  215.     endif
  216.  
  217. w = where(x lt lonmin, count)        ;Handle longitude wrap-around
  218. while count gt 0 do begin
  219.         x[w] = x[w] + 360.0
  220.         w = where(x lt lonmin, count)
  221.         endwhile
  222.  
  223.  
  224. sx = ((s1-1.)/(lonmax - lonmin))  ;Scale from lat/lon to pixels
  225. sy = ((s2-1.)/(latmax - latmin))
  226.  
  227.  
  228. ;               Now interpolate the screen image from the original.
  229. if KEYWORD_SET(bilin) THEN BEGIN
  230. ; If the image wraps around the globe, we must treat those pixels
  231. ; after the last column and before the first column specially. 
  232.  
  233.     badb = map_image_missing(image_orig, max_value, min_value)
  234.     x = (x - lonmin) * sx        ;To pixels
  235.     y = (y - latmin) * sy
  236.  
  237.     if wrap then begin
  238.     col1 = where(x gt (s1-1), count)    ;Wrap around elements
  239.     if count le 0 then goto, full_image
  240.     threecol = [Image_Orig[s1-1,*], Image_Orig[0:1,*]]
  241.     col1x = x[col1] - (s1-1)    ;Interpolate value
  242.     if badb[0] ne -1 then begin    ;Missing data value?
  243.         bad = interpolate(float(badb), x, y, miss=0)
  244.         bad[col1] = interpolate(float([badb[s1-1,*], badb[0:1,*]]), $
  245.             col1x, y[col1], miss=0)
  246.         bad = where(bad, count)
  247.         x = interpolate(Image_Orig, x, y, miss = missing) ; full image
  248.         ;Add in points that wrapped in X between s1-1 and s1.
  249.         x[col1] = interpolate(threecol, col1x, y[col1], miss = missing)
  250.         if count gt 0 then x[bad] = missing
  251.     endif else begin            ;badb
  252.         x = interpolate(Image_Orig, x, y, miss = missing) ; full image
  253.         ;Add in points that wrapped in X between s1-1 and s1.
  254.         x[col1] = interpolate(threecol, col1x, y[col1], miss = missing)
  255.     endelse
  256.     endif else begin
  257.   full_image:
  258.     if badb[0] ne -1 then begin
  259.         bad = where(interpolate(float(badb), x, y, miss = 0), count)
  260.         x = interpolate(Image_Orig, x, y, miss = missing) ;No wrap
  261.         if count gt 0 then x[bad] = missing
  262.     endif else x = interpolate(Image_Orig, x, y, miss = missing) ;No wrap
  263.     endelse
  264.  
  265. ENDIF ELSE BEGIN        ;  Scale to pixel coords & round.  
  266. ;  This is the same as: x = (x-lonmin) * sx + 0.5, but faster for arrays.
  267.     x = (x - (lonmin - .5/sx)) * sx
  268.     y = (y - (latmin - .5/sy)) * sy
  269.     if wrap then begin
  270.         bad = where(x ge s1 and x lt (2*s1), count)
  271.         if count gt 0 then x[bad]=x[bad]-s1
  272.         endif
  273.     bad = where(x lt 0 or x gt s1 or y lt 0 or y gt s2, count)
  274.     if count gt 0 then begin
  275.         x[bad] = 0
  276.         y[bad] = 0
  277.         x = Image_Orig[x, y]
  278.         badb = map_image_missing(x, max_value, min_value)
  279.         if badb[0] ne -1 then begin
  280.         worse = where(badb, count)
  281.         if count gt 0 then x[worse] = missing
  282.         endif        ;Max value
  283.         x[bad] = missing
  284.     endif else begin
  285.         x = Image_Orig[x,y]
  286.         badb = map_image_missing(x, max_value, min_value)
  287.         if badb[0] ne -1 then begin
  288.         bad = where(badb, count)
  289.         if count gt 0 then x[bad] = missing
  290.         endif        ;Max value
  291.     endelse            ;Count
  292. ENDELSE
  293.  
  294. Startx = long(screen_x[0] / scalef)
  295. Starty = long(screen_y[0] / scalef)
  296. xsize = long(nx / scalef * compress)
  297. ysize = long(ny / scalef * compress)
  298.  
  299. if compress ne 1 then $        ;Resample to screen?
  300.    x = rebin(x, nx*compress, ny*compress)  ;Interpolate screen image if necess.
  301. ;print,systime(1)-t0,' seconds'
  302.  
  303. IF scale THEN BEGIN        ;Unfudge scale factors
  304.     !x.s = !x.s / scalef
  305.     !y.s = !y.s / scalef
  306.     ENDIF    
  307. return, x
  308. end
  309.